In this lab, we will be introduced to another open source GIS software named Grass. GRASS (Geographic Resources Analysis Support System) has been under continuous development since 1982 and during the 90s a large number of federal US agencies, universities, and private companies were involved in the development. The core components of GRASS 1.0-5.0 and the management of the integration efforts into GRASS GIS releases were accomplished by the U.S. Army - Construction Engineering Research Laboratory (USA-CERL) in Champaign, Illinois. USA-CERL completed its last release of GRASS as version 4.1 in 1992, and provided five updates and patches to this release through 1995. USA-CERL also wrote the core components of the GRASS 5.0 floating point version. In 1997 and after two years of uncertainty, GRASS development was taken over by academia. Since then the international GRASS Development Team manages the source code, releases and documentation.
One of the advantages of GRASS GIS is its numerous modules for analyzing all manners of spatial and temporal data. GRASS GIS has over 500 different modules in the core distribution and over 300 addons
Included in this data are LAS files, reference maps for those files and an orthophoto of Washtenaw county. A LAS file (.las) is the standard binary format for storing airborne lidar data. The LAS dataset allows you to examine LAS files, in their native format, quickly and easily, providing detailed statistics and area coverage of the lidar data contained in the LAS files. The reference map indicates the specific location of the .las in the county. We will drape the orthophoto over our interpolated DSM.
Note: On the Mac OS, the r.in.lidar module may encounter problems. Mac user may need to use r.in.ascii in combination with libLAS instead of the provided code. 3D view may be inaccessible for Macs.
GRASS uses specific database terminology and structure (GRASS GIS Spatial Database):
A GRASS GIS Spatial Database (GRASS database) consists of a directory with specific Locations (projects) where data (data layers/maps) are stored.
Location is a directory with data related to one geographic location or a project. It is important to remember and this is a benefit of GRASS that all data within one Location has the same coordinate reference system.
A Mapset is a collection of maps within a Location, containing data related to a specific task, user or a smaller project.
New. Open the Location Wizard with the New button on the left part of the welcome screen. Select a name for the new location(this is where all your data will be stored)and create this folder.
We will establish this location’s projection-CRS (coordinate reference system). Our LiDAR data uses NAD83(HARN) / Michigan South. We will search for the correct EPSG using the Select EPSG code of spatial reference system. Use the search tool to find the correct EPSG(). After choosing a projection, a new location will be listed on the start-up screen.
Now, select the new location and mapset PERMANENT and press Start GRASS session.
In this step we will import data. This is an important step in GRASS, as it uses a specific file format only readibly in the software. Using the menu choose File - Import raster data select Common formats import In the resulting dialog box, browse to find the DEM file, change the name to DEM, and click the button Import. The imported layers will be added to GUI automatically (you can add them manually using and
.
You can also use the command line:r.import input=???.tif output=
It important to understand the computational region in Grass, which confines any analysis to a specific geographic area. The computational region can be set, subsetting a larger data extent, for quicker testing or analysis of specific regions based on administrative units. You must always set the computational region properly. These are some things keep in mind when using the computational region function:
You can check the computational region in the Console using:
g.region -p
The computational region can be set using a raster map:
g.region raster=Arb_DEM -p
Resolution can be set separately using the res parameter of the g.region module. The units are the units of the current location, in our case meters. This can be done in the Resolution tab of the g.region dialog or in the command line in the following way (using also the -p flag to print the new values):
g.region res=3 -p
The new resolution may be slightly modified in this case to fit into the extent which we are not changing. However, often we want the resolution to be the exact value we provide and we are fine with a slight modification of the extent. That’s what -a flag is for. On the other hand, if alignment with cells of a specific raster is important, align parameter can be used to request the alignment to that raster (regardless of the extent).
The following example command will use the extent from the raster named ortho, use resolution 5 meters, modify the extent to align it to this 5 meter resolution, and print the values of this new computational region settings:
g.region raster=ortho res=5 -a -p
Different functions are available through modules in GRASS (which are sometimes called tools, functions, or commands). Modules respect the following naming conventions:
| Prefix | Function | Example |
|---|---|---|
| r. | raster processing | r.mapcalc: map algebra |
| v. | vector processing | v.surf.rst: surface interpolation |
| i. | imagery processing | i.segment: image segmentation |
| r3. | 3D raster processing | r3.stats: 3D raster statistics |
| t. | temporal data processing | t.rast.aggregate: temporal aggregation |
| g. | general data management | g.remove: removes maps |
| d. | display | d.rast: display raster map |
Modules staring with v.net. deal with vector network analysis. The name of the module helps to understand its function, for example v.in.lidar starts with v so it is a vector maps function, the following name indicates module functionality importing, and finally lidar indicates that it deals with lidar point clouds. Modules and their descriptions with examples can be found the documentation. The documentation is included in the local installation and is also available online.
Here is a good tutorial showing how to manipulate and display your maps link
The fastest way to analyze a LiDAR point cloud is to use binning and create a raster map. In Grass, the r.in.lidar converts a point count (point density) into raster. As we don’t know the spatial extent of the point cloud, we cannot set the computation region ahead. Instead, we can tell the r.in.lidar module to determine the extent automatically using the -e flag. We can set the resolution (pixel size). Remember that the resolution impact the processing speed. The -r flag sets the computational region to match the output. The -o flag ignores the projection definition (the data do not have a projection, but we know that it is NAD83(HARN) / Michigan South)
r.in.lidar input=297285v.las output=Arb_297285 method=n -e -n -o resolution=10
After it has been successfully imported, we can see the spatial patterns of the LiDAR data. From the LiDAR metadata we know that the data is in feet. Grass automatically uses meters and so if you set the resolution to 3, it results in pixel size of 0.9144. We can also examine variation of the raster using the r.info function:
r.info map=Arb_297285
We can do a quick check of some of the values using the query tool in the Map Display. Examine some the pixel values
Since there are numerous points per one cell, we can use a finer resolution (notice the overwrite flag). This will result in a meter resolution map:
r.in.lidar input=297285v.las output=Arb_297285 method=n -e -n -o resolution=3 --overwrite
We can examine the distribution of the values using the histogram function. Histogram is accessible from the context menu of a layer in Layer Manager, from theMap Display`` toolbarAnalyze map``` button; or you can simply use this function in the console:
d.histogram map=Arb_297285
Now we need to do a bit of housekeeping to ensure that subsequent calculations inherit the same extent and resolution. We do this by applying the extent and resolution of our calculated image using the g.region function:
g.region raster=Arb_297285 -p
However, while this region has a lot of cells, some areas are empty due to scattering of data. To account for this absence of data in some locations, we could modify the resolution of the computational region. We use the a flag -ap to align region to resolution. For this exercise we will maintain the 1 meter resolution:
g.region res=3 -ap
A Digital Surface Model (DSM) represents the elevations of the reflective surfaces of trees, buildings, and other features elevated above the “Bare Earth”. We can use binning to obtain the DSM. Because we set the g.region, our new data will inherit the resolution and extent from the existing computational region:
r.in.lidar input=297285v.las output=binned_dsm method=max -o
To get a better understanding of different cloud values, we can compute statistics using r.report (Report raster statistics in the layer context menu):
r.report map=binned_dsm units=c
This table indicate that there are a few really high numbers (outliers at around 970 m). To better visualize this we change the color classes based on equalizing the histogram (-e flag). This allows us to see the contrast in the rest of the map (using the viridis color table):
r.colors map=binned_dsm color=elevation -e
We can also identify possible low outliers using the minimum function (method=min, we already used method=max for DSM):
r.in.lidar input=297285v.las output=minimum -o method=min
If we compare this with the range and mean of the low values, it will give an idea of how best to confine the range to eliminate possible outliers.
r.report map=minimum units=c
Again, to see the data, we can use histogram to equalized color table:
r.colors map=minimum color=elevation -e
The zrange parameter filters out points which don’t fall into the specified range of Z values. We can use these outliers to limit the range of values (748-970 m seems to be a safe bet). Use the zrange parameter to filter out any possible outliers (don’t be confused with the intensity_range parameter) when computing a new DSM:
r.in.lidar input=input=297285v.las output=binned_dsm_limited method=max -o zrange=748,970
Now we will interpolate a digital surface model (DSM) and for that we can increase the resolution to obtain as much detail as possible (we could use 1, i.e. res=1 foot, for high detail or 3, i.e. res=3 feet -a, to increase the computational speed). I will use 3 feet, but feel free to play with the resolution:
g.region raster=Arb_297285 res=3 -ap
Before interpolating, let’s confirm that the spatial distribution of the points allows us to safely interpolate. We need to use the same filters as we will use for the DSM itself:
r.in.lidar input=input=297285v.las output=count_dsm method=n return_filter=first -o zrange=748,970
Here we evaluate if there are possible outliers:
r.report map=count_dsm units=h,c,p
Then check the spatial distribution. For that we will need to use histogram equalized color table (legend can be limited just to a specific range of values: d.legend raster=count_interpolation range=0,5).
r.colors map=count_dsm color=viridis -e
First we will import the LAS file as vector points. We will only use the first return points, and limit the import vertically to avoid using outliers found in the previous step. Before running it, uncheck Add created map(s) into layer tree in the v.in.lidar dialog if you are using GUI.
v.in.lidar -t -b -o input=C:\Users\dbvanber\Box\Labs\Adv_Week_2\LiDAR_data\297285v.las output=first_returns zrange=748,970
Then we interpolate (if you have problems with the calculation you may need to set the extent using g.region raster=dsm):
v.surf.rst input=first_returns elevation=dsm tension=25 smooth=1 npmin=80
LiDAR can also be used to identify different terrain features. In this exercise we will identify different physical features of the landscape.
If you haven’t already, let’s set the computational region based on an existing raster map:
g.region raster=dsm -p
We will now import points based on the ground points only (the -t flag disables creation of attribute table and the -b flag disables building of topology; uncheck Add created map(s) into layer tree):
v.in.lidar -bt input=C:\Users\dbvanber\Box\Labs\Adv_Week_2\LiDAR_data\297285v.las output=points_ground class_filter=2 -o zrange=748,970
This figure indicates the different classes of the points in the LiDAR point cloud:
| Classification value | Meaning |
|---|---|
| 0 | Never classified |
| 1 | Unassigned |
| 2 | Ground |
| 3 | Low Vegetation |
| 4 | Medium Vegetation |
| 5 | High Vegetation |
| 6 | Building |
| 7 | Low Point |
| 8 | Reserved* |
| 9 | Water |
| 10 | Rail |
| 11 | Road Surface |
| 12 | Reserved* |
| 13 | Wire - Guard (Shield) |
| 14 | Wire - Conductor (Phase) |
| 15 | Transmission Tower |
| 16 | Wire-Structure Connector (Insulator) |
| 17 | Bridge Deck |
| 18 | High Noise |
| 19-63 | Reserved |
| 64-255 | User Definable |
Now we will interpolate the ground surface using regularized spline with tension (implemented in v.surf.rst) and at the same time we also derive slope, aspect and curvatures (the following code is one long line):
v.surf.rst input=points_ground tension=25 smooth=1 npmin=100 elevation=terrain slope=slope aspect=aspect pcurvature=profile_curvature tcurvature=tangential_curvature mcurvatur=mean_curvature
When we examine the results, especially the curvatures show a pattern which may be caused by some problems with the point cloud collection. We decrease the tension which will cause the surface to consider fewer points and increase the smoothing between points. Since the raster map called range already exists, use the overwrite flag, i.e. –overwrite in the command line or checkbox in GUI to replace the existing raster by the new one. We use the –overwrite flag shortened to just –o.
v.surf.rst input=points_ground tension=20 smooth=5 npmin=100 elevation=terrain slope=slope aspect=aspect pcurvature=profile_curvature tcurvature=tangential_curvature mcurvatur=mean_curvature --o
When we are satisfied with the result(feel free to play with the tension and smooth variables), we can compute shaded relief:
r.relief input=terrain output=relief
Now combine the shaded relief raster and the elevation raster to create a shaded relief map. Shaded relief, or hill-shading, shows the shape of the terrain in a realistic fashion by depicting how the three-dimensional surface would be illuminated from a point light source. We can calculate relief in several ways. In the GUI, we produce shaded-relief by combining the layers and changing opacity of layer or the other. We can also do this programmatically using the r.shade module, which combines the two rasters and creates a new one. We can also do this on the fly without creating any raster using the d.shade module. The module can be used from the GUI through the toolbar or using the Console tab:
d.shade shade=relief color=terrain
Alternatively, we can calculate the skyview (how it looks from an airplane) using r.skyview:
r.skyview input=terrain output=skyview ndir=8 colorized_output=terrain_skyview
Let’s combine the terrain and skyview on the fly:
d.shade shade=skyview color=terrain
We can also create a shaded relief map that highlights terrain features which wouldn’t be visible using standard shading technique. The module show relief from various directions and combines them into RGB composition:
r.shaded.pca input=terrain output=pca_shade
Local relief model. Local relief models enhance the visibility of small-scale surface features by removing large-scale landforms from the DEM:
r.local.relief input=terrain output=lrm shaded_output=shaded_lrm
The r.terrain.texture module can be used to calculate different terrain characteristics. See the the help file here. For example, we can calculate pit density:
elevation=terrain thres=0 pitdensity=pit_density
Finally, we can use the r.geomorphon module, which automatically detection of landforms (using 50 m search window):
-m elevation=terrain forms=landforms search=50
We can also isolate vegetation using LiDAR. To do this we will use the zrange parameter to filters out points which don’t fall into the specified range of Z values. Although for many vegetation-related application, a coarse resolution is necessary/appropriate due to the number of points needed to interpolate the map, we will use 1 meter:
g.region raster=dsm res=3
Let’s start by exploring the range of all the points by estimating heights for each cell. This allows us to identify tall objects infered from the difference between points at a single location:
r.in.lidar input=C:\Users\dbvanber\Box\Labs\Adv_Week_2\LiDAR_data\297285v.las output=range method=range -o zrange=748,970
To isolate and visualize tall objects, we can use r.mapcalc. This is a raster calculator in Grass. All expressions are noted using quotation marks. Here we are using a logical operator if to identify in the range map locations that are greater than 2 meters (6 feet). The cells or raster with a range greater than 2 meters are given 1, and those that are less than 2 meters are given as null:
r.mapcalc "vegetation_by_range = if(range > 6, 1, null())"
We can change the color table for a raster map from layer context menu using Set color table interactively.
Now let’s calculate the tallest objects in relation to the actual terrain. We do this by using the max points and using the base_raster we already calculated. The function assumes that the base map is zero and we set the range to 100 meters as it is unlikely that there are such tall objects
r.in.lidar -d input=C:\Users\dbvanber\Box\Labs\Adv_Week_2\LiDAR_data\297285v.las output=height_above_ground method=max base_raster=terrain zrange=0,100 -o
Now let’s redo the calculation to identify locations with objects that are higher than 2 m. If we wanted to preserve the elevation in the output, we could use if(height_above_ground > 6, height_above_ground, null()) in the expression, but we want just the areas, so we use constant (1).
r.mapcalc "above_2m = if(height_above_ground > 6, 1, null())"
Often there are speckles where there is no data due to issues with the point cloud interpolation. To reduce these speckles, we can use the r.grow, which assign data around existing non-null cells/rasters:
r.grow input=above_2m output=vegetation_grow
This is likely a good approximation of locations with vegetation and the occasional building (which we can filter out later). We calculate the clump of individual cells into patches using r.clump. This helps identify contiguous area of forests and filter out small insignificant vegetation:
r.clump input=vegetation_grow output=vegetation_clump
Some of the patches are very small. Using r.area we remove all patches smaller than the given threshold:
r.area input=vegetation_clump output=vegetation_by_height lesser=100
Now convert these areas to vector:
r.to.vect -s input=vegetation_by_height output=vegetation_by_height type=area
So far we were using elevation of the points, now we will use intensity. Intensity is a measure, collected for every point, of the return strength of the laser pulse that generated the point. It is based, on the reflectivity of the object struck by the laser pulse intensity is used by r.in.lidar when -j (or -i) flag is provided:
r.in.lidar input=C:\Users\dbvanber\Box\Labs\Adv_Week_2\LiDAR_data\297285v.las output=intensity zrange=748,970 -j -o
Given this high resolution, there are some cells without any points, so we use r.fill.gaps to fill these cells (r.fill.gaps also smooths the raster as part of the gap-filling process).
r.fill.gaps input=intensity output=intensity_filled uncertainty=uncertainty distance=3 mode=wmean power=2.0 cells=8
Let’s apply a grey scale color table as it is more appropriate for intensity:
r.colors map=intensity_filled color=grey
There are some areas with very high intensity. We can equalize the color table using the histogram to make it easier to visually explore the data
r.colors -e map=intensity_filled color=grey
Let’s use r.geomorphon again, but now with DSM to identify specific features that can be inferred from the intensity layer. The following settings shows building roof shapes:
r.geomorphon elevation=dsm forms=dsm_forms search=7 skip=4
With different settings, especially a higher flatness threshold (flat parameter), we can isolate forested areas. Using the geomorphon function vegetation appears as footslopes, slopes, and shoulders. Individual trees are represented as shoulders.
r.geomorphon elevation=dsm forms=dsm_forms search=12 skip=8 flat=10 --o
Lowering the skip parameter decreases generalization and brings in summits which often represents tree tops:
r.geomorphon elevation=dsm forms=dsm_forms search=12 skip=2 flat=10 --o
Now we extract summits (trees) using r.mapcalc raster algebra expression if(dsm_forms==2, 1, null()).
r.mapcalc "trees = if(dsm_forms==2, 1, null())"
How would you get a vectors of tree locations?
We can explore our study area in the 3D view:
We might also want to mask the data to a specific location using r.mask(this is clip in many gis platforms). You can download .las file for Ann Arbor here() for doing your own terrain models. Ensure that you choose the correct computational surface.
You can execute Python code in the Simple Python Editor (accessible from the toolbar or the Python tab in the Layer Manager). Another option is to use your favorite text editor and then run the script in GRASS GIS using the main menu File → Launch script.
We can use the Simple Python Editor to run commands. When you opening the Editor, you find a short code snippet that imports the GRASS GIS Python Scripting Library:
import grass.script as gscript
You can run Grass function using the gscript.run_command. For example, if we want to call the g.region function we run this code:
gscript.run_command('g.region', flags='p')
In this example, we set the computational extent and resolution to the raster layer ortho which would be done using g.region raster=ortho in the command line. To use the run_command to set the computational region, replace the previous g.region command with the following line: